Assessment

How did antibiotic prescription trends vary pre-, mid-, and post-lockdown in different regions of Scotland?

#How did the prescription of the top 10 most prescribed antibiotics in 2019 change due to lockdown

#How did covid-19 pandemic and lockdown affect the prescription trends of the most common antibiotics in Scotland. (Would it be better to focus on the main lockdown (May 2020) or also compare between lockdowns?)

First I will load the appropriate libraries. Then, load, clean and filter the data.

First, I loaded the datasets and created a function to concisely use the same functions on different datasets. I learned how to do this from the website https://www.dataquest.io/blog/write-functions-in-r/. In this function I filtered for antibiotics by only looking at prescriptions with the BNFItemCode starting with 0501. Then I filtered out NAs, cleaned up the names and calculated the number of prescriptions. I then merged the datasets

may2020 <- read_csv(here("data","pitc202005.csv"))
may2024 <- read_csv(here("data","pitc202405.csv"))
may2019 <- read_csv(here("data","pitc201905.csv"))

may2019 <- may2019 %>% 
  rename("HBT" = HBT2014)

process_antibiotics <- function(data) {
  data %>%
    filter(str_detect(BNFItemCode, "^0501"), !is.na(BNFItemDescription)) %>% 
    mutate(BNFItemDescription = str_replace(BNFItemDescription, "_", " "),
           Antibiotic = str_extract(BNFItemDescription, "^[^ ]+")) #Extracts the first word, which is the antibiotic name
}
#I filtered for prescriptions with an item code beginning with 0501 because prescriptions with that item code are antibiotics (Open Prescribing, 2019).

may2019_processed <- process_antibiotics(may2019)
may2020_processed <- process_antibiotics(may2020)
may2024_processed <- process_antibiotics(may2024)

#I split this into 2 functions as later I want to group by HBT
antibiotics_by_name <- function(data2)  {
  data2 %>% 
    group_by(Antibiotic) %>%
    summarise(quantity_sum = sum(PaidQuantity), .groups = "drop") %>%
    arrange(desc(quantity_sum))
}

#I filtered for prescriptions with an item code beginning with 0501 because prescriptions with that item code are antibiotics (Open Prescribing, 2019).

may2019_processed <- antibiotics_by_name(may2019_processed)
may2020_processed <- antibiotics_by_name(may2020_processed)
may2024_processed <- antibiotics_by_name(may2024_processed)

I then prepared the data for joining by renaming Fluclox to Flucloxacillin. I then joined the data, and tidied it

may2020_processed <- may2020_processed %>% 
  mutate(Antibiotic = ifelse(Antibiotic == "FLUCLOX", "FLUCLOXACILLIN", Antibiotic))

may2019_processed <- may2019_processed %>% 
  mutate(Antibiotic = ifelse(Antibiotic == "FLUCLOX", "FLUCLOXACILLIN", Antibiotic))

antibiotics1920 <- full_join(may2019_processed, may2020_processed, by = "Antibiotic")
antibiotics192024 <- full_join(antibiotics1920, may2024_processed, by = "Antibiotic")

antibiotics192024 <- antibiotics192024 %>% 
  rename(
   "May 2019" = quantity_sum.x,
   "May 2020" = quantity_sum.y,
   "May 2024" = quantity_sum) %>% 
  slice(1:10)

antibiotics192024 <- antibiotics192024 %>% 
  mutate(Antibiotic = str_to_sentence(Antibiotic)) #I changed all the antibiotic names to only have the first letter capitalised. I learnt this from the R help page

I created a stacked bar graph to visually compare prescriptions

pivoted_antibiotics <- antibiotics192024 %>% 
  pivot_longer(
    cols = c("May 2019", "May 2020", "May 2024"), 
    names_to = "Year", 
    values_to = "Prescription_Number")

options(scipen=999)

antibiotics_graph <- pivoted_antibiotics %>% 
  ggplot(aes(x = Prescription_Number, y = Year, fill = Antibiotic)) +
  geom_bar(stat = "identity") +
  labs(title = "Prescriptions Sold for Each Antibiotic by Year") +
  scale_fill_brewer(palette = "Paired") +
  theme_minimal()

antibiotics_graph <- antibiotics_graph %>% 
  ggplotly(tooltip = c("Prescription_Number", "Antibiotic"))
antibiotics_graph

plotly - interactive graphs in R

This data visulaisation shows that the overall number of prescirptions sold decreased during lockdown and was slightly lower before lockdown than after. It also shows that the top 10 most prescribed antibiotics were not the same in every year.

I then created a gt() table to show the number of prescriptions and calculate the total prescriptions per year

all_antibiotics_table <- antibiotics192024 %>%
  summarise(across(starts_with("May"), sum, na.rm = TRUE)) %>%  
  mutate(Antibiotic = "Total") %>%
  bind_rows(antibiotics192024, .)

all_antibiotics_table <- all_antibiotics_table %>%
  gt() %>% 
  cols_label(Antibiotic = "Antibiotic") %>% 
  cols_align(align = "center",
             columns = everything()) %>% 
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels(columns = everything())
  ) %>% 
  tab_header(
    title = "Affect of the Covid19 pandemic on the prescription of common antibiotics",
    subtitle = "Looking at the change of prescriptions for the top 10 most common antibiotics from before, during and after the longest Covid19 lockdown") %>% 
  tab_spanner(
    label = "Number of Prescriptions",
    columns = c("May 2019", "May 2020", "May 2024"))
all_antibiotics_table

Antibiotic	Number of Prescriptions
Affect of the Covid19 pandemic on the prescription of common antibiotics
Looking at the change of prescriptions for the top 10 most common antibiotics from before, during and after the longest Covid19 lockdown
Antibiotic	May 2019	May 2020	May 2024
Amoxicillin	2875868	1778491.0	3751932.0
Phenoxymethylpenicillin	1892108	968754.0	2616743.0
Flucloxacillin	1612959	1475567.0	1733244.0
Trimethoprim	782422	745904.0	696670.0
Erythromycin	661545	443204.0	582068.0
Lymecycline	623567	493844.0	538300.0
Oxytetracycline	572089	428371.0	277441.0
Doxycycline	447389	387163.0	616170.0
Co-amoxiclav	382770	349958.0	466694.0
Clarithromycin	377144	230282.1	493464.5
Total	10227861	7301538.1	11772726.5

I then created a faceted geospatial map to show which areas of Scotland were prescribed more antibiotics than other areas, and how this changed over covid

#Load data
NHS_healthboards <- st_read(here("data", "SG_NHS_HealthBoards_2019.shp"))

## Reading layer `SG_NHS_HealthBoards_2019' from data source 
##   `C:\Users\franc\OneDrive\Documents\Data Science\Assessment\B207252\data\SG_NHS_HealthBoards_2019.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 5512.998 ymin: 530250.8 xmax: 470332 ymax: 1220302
## Projected CRS: OSGB36 / British National Grid

#Prep for joining
NHS_healthboards <- NHS_healthboards %>% 
  rename("HBT" = HBCode)

#Summarise antibiotic data by HBT code to give a total per area
antibiotics_by_HBT <- function(data_frame) {
    data_frame %>% 
      group_by(HBT) %>%
      summarise(quantity_sum = sum(PaidQuantity), .groups = "drop") %>%
      arrange(desc(quantity_sum))
}

#Apply the function
HBTmay2019_processed <- antibiotics_by_HBT(may2019)
HBTmay2020_processed <- antibiotics_by_HBT(may2020)
HBTmay2024_processed <- antibiotics_by_HBT(may2024)

#Join antibiotic data
HBTantibiotics1920 <- full_join(HBTmay2019_processed, HBTmay2020_processed, by = "HBT")
HBTantibiotics192024 <- full_join(HBTantibiotics1920, HBTmay2024_processed, by = "HBT")

#Rename and remove NAs
HBTantibiotics192024 <- HBTantibiotics192024 %>% 
  rename(
   "May 2019" = quantity_sum.x,
   "May 2020" = quantity_sum.y,
   "May 2024" = quantity_sum) %>% 
  filter(!if_any(everything(), is.na))

#Pivot to make faceting easier
HBTpivoted_antibiotics <- HBTantibiotics192024 %>% 
  pivot_longer(
    cols = c("May 2019", "May 2020", "May 2024"),
    names_to = "Year", 
    values_to = "Prescription_Number")

#Join with healthboard data

Healthboards_and_antibiotics <- NHS_healthboards %>% 
  full_join(HBTpivoted_antibiotics, by = "HBT")

#Geom plot
Healthboards_and_antibiotics_plot <- Healthboards_and_antibiotics %>% 
  ggplot(aes(fill = Prescription_Number)) +
  geom_sf(color = "white", linewidth = 0.1) +
  scale_fill_viridis_c(name = "Total No. of Prescriptions",
                       labels = scales::comma) +
  labs(title = "Total Number of Antibiotic Prescriptions in Areas of Scotland from 2019 to 2024") +
  facet_wrap(~ Year) +
  annotation_scale(    
    location = "tl"
  ) +
  annotation_north_arrow(
    location = "tl",    
    pad_y = unit(0.5, "in"),    
    style = north_arrow_nautical(
      fill = c("grey40", "white"),      
      line_col = "grey20"
    ) 
  )

geog_graph <- Healthboards_and_antibiotics_plot %>% 
  ggplotly()

Assessment

Francesca

2024-11-01

How did antibiotic prescription trends vary pre-, mid-, and post-lockdown in different regions of Scotland?